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Abstract: We present the first fully automated implementation of cross-section com¬ 
putation and event generation for loop-induced processes. This work is integrated in the 
MadGraph5_aMC@NLO framework. We describe the optimisations implemented at the 
level of the matrix element evaluation, phase space integration and event generation allow¬ 
ing for the simulation of large multiplicity loop-induced processes. Along with some selected 
differential observables, we illustrate our results with a table showing inclusive cross-sections 
for all loop-induced hadronic scattering processes with up to three hnal states in the SM 
as well as for some relevant 2—7-4 processes. Many of these are computed here for the hrst 
time. 
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1 Introduction 

The first run of the Large Hadron Collider (LHC) further confirmed the Standard Model 
(SM) in a spectacular way with the discovery of what was its last missing piece: the Higgs 
boson [1-5]. This discovery testifies not only of the technical success of the accelerator 
and detectors, but also of the level of accuracy reached in the theoretical predictions in¬ 
tervening at many levels of a discovery or exclusion at a collider. One can identify three 
main types of effort at the origin of this precision reached by modern high energy physics 
simulations. First, dedicated analytical single-purposed computations keep pushing the 
boundaries of perturbative physics and are essential for realistic inclusive predictions {e.g. 
Higgs hadroproduction at N^LO [6]). Secondly, parton shower Monte-Carlo programs such 
as HerwigH —h [7], Pythias [8] and Sherpa [9] improved their formal control on the soft 
physics resummation as well as on the merging techniques with hard matrix-element predic¬ 
tions. Finally, the last decade has seen the rise of a number of automated one-loop matrix 
element computation tools, such as MadLoop [10], FeynArts [11, 12], OpenLoops [13] 
and GoSam [14]. Thanks to these progresses, one can obtain accurate normalisation and 
distributions for basically any process of interest at NLO accuracy, using flexible partonic 
hard event generators such as POWHEG [15-17], Sherpa and the framework in which 
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the present work is carried, MadGraph5_aMC@NLO [18] (referred to as MG5 aMC hence¬ 
forth) . 

The aforementioned tools are crncial for both experimentalists and model builders to 
simulate backgrounds and explore new physics signals. However, despite these efforts, event 
generation for loop-induced processes, i.e. processes without any tree-level contributions 
and starting at one-loop, is still not available in a systematic and fully automated way. This 
is to be contrasted with the fact that loop-induced amplitudes for many relevant processes 
are already available in the main public one-loop matrix element providers [10-14]. The 
aim of this work is to remedy this situation by providing an efficient and generic technique 
for the simulation of loop-induced processes which play a significant role both in the SM 
and beyond. Loop-induced processes take their name from the fact that their leading- 
order (LO) contribution comes already from loop amplitudes, so that loop integrals are 
unavoidable in this case, even for the crudest computation of the cross-section. Exceptions 
to this are certain processes where the loop featuring heavy particles can be integrated out 
to form an effective point-like vertex turning the loop topology into a tree one. However, 
this approximation is typically valid only in a limited kinematical range, so that the general 
implementation of the direct computation of loop-induced processes is desirable. 

In the SM, a prime example is the gluon fusion channel for Higgs boson production 
which dominates the inclusive production, mainly because of the large gluon luminosity 
in high energy hadron colliders. Besides this obvious case, loop-induced channels can also 
amount to a significant part of some NNLO corrections; it was for example recently shown 
in ref. [19] that the contribution of the loop-induced gg —?■ ZZ channel represents 60% of 
the full NNLO correction of hadronic Z-boson pair production. Loop-induced processes are 
also often relevant in the context of beyond the Standard Model (BSM) models where the 
associated loop suppression factor can help to evade current experimental bounds. 

The attitude towards the computation of loop-induced processes has been so far that 
of developing single-purpose specific codes (see references in sect. 3). Given current loop 
technology, it is desirable to adopt a solution both generic in the process studied and flexible 
in its use. We present here the implementation of this solution within the public frame¬ 
work MGSaMG. Loop matrix-elements are computed with MadLoop, using a combination 
of the Ossola-Papadopoulos-Pittau (OPP) [20, 21] reduction method as implemented in 
GutTools [22] and various Tensor Integral Reduction tools (TIR) such as PJFry [23], 
Golem95 [24, 25] and the in-house implementation IREGI. The integration over phase- 
space and event generation is performed by MadEvent [26], the multi-purpose phase space 
integrator used for LO calculations in MGSaMG. 

The paper is organised as follows. Sect. 2 describes the various improvements brought 
to MadLoop and MadEvent to cope with loop-induced processes, for which the absence 
of underlying tree topologies renders inapplicable many optimisations. Sect. 3 presents 
a comprehensive list of cross-sections for SM loop-induced processes obtained within our 
implementation. In sect. 4, we focus on loop-induced Higgs production to illustrate how the 
various simulation features of MGSaMG apply in this case. A systematic comparison of our 
results with those obtained using the Effective Field Theory (EFT) approach shows where 
and to which extent the latter is a good approximation. Sect. 5 is devoted to the validation of 
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our implementation within a BSM context against the independent computation of ref. [27] 
for loop-induced Z-Higgs associated production. We also show results for the loop-induced 
charged Higgs pair production. We summarise our work in sect. 6. 


2 Technique 


2.1 Loop computation, TIR implementation and colour decomposition 

The main challenge for loop-induced processes integration is the running speed of the related 
loop-matrix element computation. In NLO computations, the bulk of the contribution 
comes from real-emission and Born topologies, so that only limited statistics is necessary 
for the evaluation of the virtual contribution (see sec. 2.4.3 of ref. [18]); it is then very 
often not the limiting factor in terms of computational load. The situation is radically 
different for the direct integration of loop-induced processes where the computation time is 
proportional to the execution speed of the loop matrix element. 

To understand the characteristics of the various techniques available in MadLoop for 
the computation of loop-induced squared matrix elements, it is appropriate to start from 
its generic expression: 
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where designates the loop-induced amplitude. A/, are colour structures, are 

loop integrand numerators and Di,ij are d-dimensional propagator denominators of the form 
— mf,. The symbols X^coiour denote the sum over all helicity, colour 

configurations and loop subamplitudes factoring a single colour factord The integers i/, 
L and are the total number of helicity configurations, loop subamplitudes and loop 
propagator denominators in subamplitude It respectively. Notice that no UV counterterms 
are necessary in this case given that loop-induced processes are finite. Conversely, R2 
counterterms [28-30], originating from the 4-dimensional nature of the reduction methods 
applied in MadLoop, must be included.^ The \ term of eq. 2.1 is the new element 

specific to loop-induced matrix element computations, and we now turn to detailing its 
implementation in MadLoop. 

The key characteristic of the quantity is that it is not linearly dependent on 

the loop amplitudes. This has dramatic implications on the type of optimisations applicable 
to its computation. When replacing the loop integral by the formal reduction operator Red[] 
(symbolically denoting the application of any one of the available reduction tools interfaced 

^ These subamplitudes are in one-to-one correspondence with the amplitudes of the constituting Feynman 
diagrams, except for those involving vertices featuring mlultiple colour factors {e.g. four-gluon vertex). 

^Counterterms are selected and constructed with the same algorithm as for NLO virtual matrix element 
generation, i.e. directly starting from the loop Feynman diagrams considered, hence guaranteeing the 
consistency of the computation. 
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to MadLoop), one finds: 
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Contrary to NLO computations, where loop amplitudes interfere with tree-level ones, it is 
clear that the Red[] operator cannot be pulled out of the sum over helicity conhgurations so 
as to apply at the squared matrix element level instead (see the transition from eq. 2.75 to 
eq. 2.76 of ref [18]). The expanded sum contains terms and this quadratic 

scaling with the number of diagrams is problematic when compared to the apparent^ linear 
scaling of the total loop computation time. As a result, the computation time becomes 
dominated by the squaring operation for loop-induced processes with as few as a thou¬ 
sand diagrams. To circumvent this scaling, we consider here the same solution as the one 
adopted for tree-level computations, namely colour decomposition. Indeed, the origin of 
this problem can be traced back to the size of the colour basis, built out of a total of L basis 
vectors A; which are not linearly independent. The solution consists then in projecting the 
colour factors A; onto a colour-flow basis [31] built out of K colour-flow basis vectors Ki 
(chains of Kronecker delta structures with indices in the (anti-)fundamental representation 
of SU{3)c)- The growth of K with the multiplicity is power-like, hence guaranteeing that 
the computational cost of the colour algebra involved in the amplitude squaring operation 
remains negligible with respect to that of the loop amplitude computation. For example, 
the process gg —)• hggg has 3330 subamplitudes but only 24 different colour-flows. Each 
loop colour factor is then projected onto the colour-flow basis as follows: 
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Notice that in practice, the projection matrix is sparse so that most of the projection 
coefficients ai^i are zero and the sum involves a few terms only. The corresponding colour 
matrix is 

^ k,k] = Kij. (2.3) 

colour 

Once expressed in the colour-flow basis, the loop-induced squared amplitude reads 

= EEE (2-t) 

h=li=l j=l 

®This statement must be taken with care since both the construction and reduction of loops increases 
in complexity with the multiplicity. However, internal recycling of currents effectively implements recursive 
relations, mitigating this increased complexity and yielding a scaling approximately linear in the number 
of diagrams, at least for processes with up to four Hnal state legs [13]. 
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where Ji^h is defined as follows: 
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and corresponds to the partial colour amplitudes built from sums of the Lorentz part of the 
loop amplitudes weighted by projection coefficients. This change of colour basis is not only 
computationally advantageous, but also more physical since the partial colour amplitudes 
are gauge invariant. It also automatically solves the problem of assigning a colour flow 
to partonic events, which amounts to specifying the colour dipole pairs in the starting 
conditions of parton shower Monte-Carlo programs and thus has an important impact on 
the radiation pattern. 

This colour projection is performed automatically by MG5aMC, using the same colour 
algebra module employed for tree-level matrix element generation. However, contrary to 
the tree-level case, the projection coefficients ai^i can be of different orders in the colour 
expansion in because of the loop colour trace. This effect is accounted for when assigning 
colours to the generated events in which case only the leading term in the colour expansion 
is kept. Also, thanks to the flexibility of the colour module performing the algebra, arbitrary 
colour structures can be supported as well as the definition of any other basis, were that be 
necessary. The computation of the partial colour amplitudes has been made available also 
for the computation of the standard loop and Born interference term appearing in NLO 
computations (see appendix A. 2). 

The inability to perform loop reduction at the squared matrix element level leads to a 
crucial difference between the TIR and OPP reduction methods. To understand why this 
is so, we detail the expression taken by the Red[] operation in both cases. The integrand 
numerator is decomposed in a polynomial in the loop momentum i as follows: 

( 2 . 6 ) 
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where ^ ^ are referred to as the polynomial coefficients. In the case of OPP reduction, 

this decomposition only serves the purpose of improving the efficiency of the computation 
since the polynomial coefficients can be recycled for all the different values of the loop 
momentum for which N{t)i^h must be evaluated. In the case of TIR, this decomposition is 
essential since the object reduced are the tensor integrals ■. 
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We can now write the more precise form taken by OPP and TIR reduction: 
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It is now manifest that the output of the OPP reduction depends on both the loop 
and helicity considered while the TIR output only depends on the loop considered [i.e. the 
tensor integrals ^ only carry a dependence on the index I, not h). For this reason, the 
number of independent OPP reductions performed per kinematic configuration is necessar¬ 
ily LxH, that is the number of loop amplitudes times the number of helicity configurations. 
On the other hand, the number of independent^ TIR per phase-space points is only pro¬ 
portional to L since tensor integrals can be recycled across helicity configurations. For this 
reason, the evaluation of loop-induced matrix elements summed over helicity configurations 
is typically faster using TIR. Conversely, for the computation of a single helicity config¬ 
uration, OPP type of reduction is preferable since the implementations of TIR currently 
available in MadLoop, namely PJFry, Golem95® and the in-house implementation IREGI, 
are slower for an equal number of calls to the Red[] operator (see appendix B for quanti¬ 
tative results on benchmark processes). Note that this statement is highly dependent on 
the TIR implementation considered, and comparisons presented in ref. [13] suggest that it 
might not hold true for the private TIR implementation in GOLLIER [32, 33]. However, 
due to the different scaling of the complexity of TIR and OPP reduction with the rank of 
the loop integral, it is expected that OPP is always faster for larger ranks (typically r > 6). 

Given the above, it is unclear whether event generation is more efficient using TIR 
and an explicit sum over helicity configurations or using OPP reduction and a Monte-Carlo 
(MC) sampling over them. The outcome mostly depends on the quality of the helicity 
discrete importance sampling (i.e. the magnitude of its mild dependence with kinematics), 
the runtime speed of TIR implementations and the relative importance of the computational 
cost of the determination of the polynomial coefficients (which scales like Lx H). Our tests 
show that generating events with an MC over helicity configurations, using an independent 
importance sampling for each integration channel, yields better timings and we use it as 
our default (see sect. 2.2). 

We now turn to the description of two additional minor improvements on MadLoop. 
First, many loop-induced processes occur only via closed fermion loops where often the 
different massless flavours bring the exact same contribution and can therefore be recycled. 
Identical diagrams (couplings, propagator spins, masses and widths are compared) are de¬ 
tected at generation level and then traded for a multiplicative factor affecting one chosen 
representative diagram.® For example, this leads to an improvement of a factor two for the 
production of electroweak bosons via gluon fusion. Secondly, because of Furry theorem, 
Feynman loop diagrams with an odd number of photon external legs can be exactly zero 
and their presence slows down the integration because MadEvent must probe the corre¬ 
sponding channels many times before deciding that it can be safely discarded. To avoid this, 
MadLoop has been modified so as to detect and remove such Furry loops at the diagram 

^They are not completely independent, since most TIR implementations internally cache some interme¬ 
diate reduction results and scalar integral computations so that they can be re-used across loops sharing 
some of their reduced topologies. 

®The current version of Golem 95 does not currently allow tensorial coefficients to be recycled. 

®This means that the loop diagrams removed this way no longer show in the list of diagrams drawn by 

MadLoop. 
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generation step (with a notification to the user). 

The loop-induced matrix element codes generated by MadLoop have been validated in 
several ways. First, we compared numerical results for specihc kinematic conhgurations of 
many processes against MadLoop4 [10], whose implementation is completely independent 
of MG5aMC, and against MadLoopS in a mode that does not use the polynomial decom¬ 
position of the integrand numerator (see appendix A. 2). The list of processes cross-checked 
includes, among others, dd —?■ hggg and all gluon fusion processes with up to three identical 
neutral massive bosons in the hnal states. To facilitate future comparisons, we report in 
appendix B the numerical result for a chosen phase-space point of the process gg —)• hhgg. 
Secondly, at the integrated level, we compared our prediction for the partial decay width 
z —)• ggg to the result reported in ref. [34] . ' Additional comparisons with the Higgs Effective 
Theory are presented in sect. 4. Finally, we tested Lorentz invariance, crossing symmetries 
and gauge invariance using Ward identities with the standard checks implemented via the 
command ‘check’ of MGSaMC interface. 

2.2 Phase-space integration and event generation 

The phase-space integration and event generation is based on the MadEvent [26] algorithm 
that we improved in the context of this work. At its core lies the diagram enhancement 
method which separates the integration into a sum of integrals whose singularity structure is 
dictated by a single Feynman diagram topology. Each of these individual integrals, referred 
to as integration channel, is then integrated using an appropriate phase-space parametri- 
sation undoing its underlying structure. For this reason, using MadEvent phase-space 
mapping algorithms requires to build tree-level topologies from the contributing loop dia¬ 
grams. This is achieved by contracting the loop to a single vertex point. It is important 
to stress here that this mapping to tree topologies is only used to setup the phase-space 
parametrisation, and at no point to build any sort of numerical estimate of the correspond¬ 
ing loop-induced amplitude (the full exact loop-induced matrix element is used throughout 
our implementation). 

A MadEvent run involves two steps. The hrst one is referred to as the survey and 
consists in computing the cross-section for each integration channel down to a given accu¬ 
racy of 5%. Using the information on relative cross-sections and efficiencies provided by 
the survey, MadEvent proceeds with a second step referred to as the refine where event 
generation takes place. The time necessary to run the rehne step is approximatively linearly 
proportional to the requested number of unweighted events whilst the one of the survey is 
independent of this number. 

We have implemented a series of improvements to MadEvent, the most important of 
which being the implementation of a dynamical importance sampling for the Monte-Carlo 
over helicity conhgurations (he. the frequency of probing a particular helicity conhguration 
is dynamically adjusted to its relative contribution). Notice however that, as it is the case 
for standard adaptive Monte-Carlo [35], we do not account for any correlation between the 

^The result reported in item g.7 of table 8 differs from the one of ref. [34] where the value of the top 
mass, as and a~^ are chosen to be 173 GeV, 0.134 and 128.0 respectively. 
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Helicity sum 

Monte-Carlo 

Exact 

Loop Reduction 

CutTools 

CutTools 

TIR 


Survey 


pp hj 

13m (125k) 

32m (260k) 

9m (260k) 

pp hjj 

2d4h (1.2M) 

16dl0h (5.4M) 

9dl3h (5.4M)* 

gg zz 

lh06m (34k) 

12h50m (255k) 

lh44m (255k) 

gg zhg 

llhl3m (110k) 

ld8h (516k) 

ld4h (516k)* 

Refine 

pp —)• hj 

lh43m (385k) 

23m (431k) 

6m (431k) 

pp hjj 

7dl7h (2.18M) 

75dlh (20.6M) 

51dl9h (20.6M)* 

gg zz 

7h20m (407k) 

4dl3h (4.55M) 

23h07m (5.78M) 

gg zhg 

23h03m (277k) 

2d22h (1.13M) 

3dl4h (1.4M)* 


Table 1. Cumulated CPU-hours necessary for the generation of 10k unweighted events for various 
loop-induced processes. This corresponds to a Monte-Carlo accuracy on the cross-section typically 
better than a percent. The numbers in parenthesis specify the number of polarised (even when 
summing exactly over helicity configurations) phase-space points for which the matrix elements 
needed to be evaluated. The column ‘TIR’ only reports the fastest timing between using the 
PJFRY and IREGI implementation of TIR and we suffixed the timing by a star (*) when the 
latter is faster. For each process, we underlined the fastest of all approaches by reporting it in bold 
font. 


kinematic variables of integration and the sampling distribution of helicity configurations. 
The gain obtained thanks to this new operational mode is made explicit in table 1. In that 
table and for a small set of loop-induced processes, we present the cumulated CPU-hours^ 
necessary for the survey step as well as for the refine to 10k unweighted events. We also 
report the total number of phase-space points for which the matrix-elements needed to be 
evaluated during the integration.® We warn the reader that the quantitative results reported 
here must not be interpreted too literally given that the speed of the node assigned by the 
cluster can vary from one run to the other. We compare timings for an integration using a 
Monte-Carlo sampling with an explicit sum over helicity configurations. In the latter case, 
we further compare three different reduction methods; OPP as implemented in CutTools 
and TIR as implemented in PJFry or IREGI. 

We find that when summing exactly over helicity configurations for each phase-space 

®This is equivalent to how long the integration would have taken if it was run sequentially on a single 
CPU. 

^Notice that for the processes pp —>■ hj and pp —>■ hjj, the repartition of phase-space points probed 
across different partonic subprocesses is not the same for the survey and refine steps or between the two 
different integration techniques. This explains for example the difference of timing per phase-space point 
between the survey and refine steps. 
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points, TIR is advantageous thanks to the recycling of tensorial coefficients across these 
configurations (see sect. 2.1). However, in the case of processes of larger multiplicity (pp —>• 
hjj) or with many helicity configurations {gg —?• zz), this is not sufficient to overcome 
the gain obtained from sampling over helicity, which considerably reduces the number of 
polarised phase-space points that needs to be probed to reach a given accuracy. 

Given that the typical time for the evaluation of a 2 —t- 4 loop-induced matrix element 
for a phase-space point can already reach several seconds (see table 14 in appendix B), it 
is of paramount importance to be able to scale MadEvent parallelisation independently 
of the number of integration channels. In this way, even though the cumulated sequen¬ 
tial CPU-hours necessary for the computation remains constant, the user time spent for 
the integration of a process can be reduced proportionally to the number of cores available. 
Although Monte-Carlo integration methods are intrinsically trivially parallel, adaptive tech¬ 
niques, where the phase-space probing distribution is improved over time, break the mutual 
independence of two successive iterations. The original approach of MadEvent implements 
the simplest parallelisation in this context, which consists in running each channel of in¬ 
tegration as an independent separate job. However, this is not sufficient for loop-induced 
processes which exhibit a small number of time-consuming channels. To push parallelisa¬ 
tion further, we have therefore modified the steering of MadEvent so as to submit multiple 
jobs on the cluster for each iteration. In between iterations, the phase-space sampling grids 
and cross-section results from each job are combined together to build the input of the next 
iteration which is again split into a series of new independent jobs. 

Finally, we have also improved event generation. In MadEvent, the unweighting op¬ 
eration follows each iteration and if it did not generate enough unweighted events then 
a new iteration is submitted -setup to probe twice as many phase-space points- and the 
events already generated are simply discarded. We ameliorated the unweighting step by 
allowing to combine events generated in any of the previous iterations. One caveat lies in 
the fact that if the MC sampling grids improve significantly from one iteration to the other, 
this procedure does not take advantage of the resulting increase in unweighting efficiency. 
We remedy this problem by estimating the remaining computing time for each of the two 
strategies (discarding or not events from previous iterations) and choosing the most efficient 
one. 

Using a still private implementation of an interface between MadLoop and SHERPA, 
we carried a detailed comparison at the integrated level for the loop-induced process 
gg —)• ZZ, including both the triangle (featuring an s-channel Higgs) and box topologies. 
In particular, we compared to the per-mil accuracy the cross-sections for each individual 
squared topology as well as for their interference only. We found perfect agreement between 
MadEvent and SHERPA. We also checked that integrating all contributions together yields 
consistent results. Notice that in both computations of this check, the loop-induced ma¬ 
trix element implementation is identical, being the one provided by MadLoop, and as 
such this validation targets only MadEvent integration procedure. Additional validations 
of the phase-space integration are presented in sect. 4 and 5. Finally, collaborators have 
compared our predictions with results for the processes pp —)• VV obtained completely 
independently [36]. 
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3 Results in the Standard Model 


In this section, we present results for inclusive cross-sections of loop-induced processes 
within the Standard Model, whose parameters are set to the values listed in table 2. 


Parameter value 


Parameter value 


as(m|) 

fJ-R = IJ‘F 
mt = yt 
Gf 
mz 

mw 

mu 

yCKM 
mr± = yr± 


I set by PDF if present (0.13355) 
I 0.118 otherwise 

p. (see text for def.) 

173.0 

1.16639e-05 

91.188 

Mz /1 /i 47r a 

V2Y V V2GPM^ 

125.0 

Sij 

1.777 


nif 

4 

mb = yb 

4.7 

Tt 

0 

a~^ 

132.507 

Tz 

2.4414 

Tw 

2.0476 

Th 

0.00638 

me± = m^± 

0.0 

r.± 

0.0 


Table 2. Standard Model parameters used for obtaining the results presented in tables 4-8. 
Dimensionful parameters are given in GeV. 


We considered all hadronic loop-induced processes in the SM up to three particles in 
the final states, featuring one, two or three heavy bosons and/or photons as presented in 
tables 4, 5 and 6 respectively. The seemingly missing processes are all zero, either because of 
Furry’s theorem (such as in gg —?• ha for instance, where all fermionic loops with clockwise 
and anti-clockwise flow cancel pair-wise) or Landau-Yang’s theorem [37, 38] that forbids a 
massive vector particle to decay into two identical massless ones, such as in gg —?■ z. In 
addition, we present in table 7 cross-sections obtained for a selected list of processes with 
four external final states as well as loop-induced processes with non-hadronic initial states. 

For each process, we generated a sample of 10k unweighted events, yielding a Monte- 
Carlo accuracy on the inclusive cross-section of at least 1% but typically better. The central 
factorisation and renormalisation scales are set dynamically to half the sum of all final state 
transverse energies in the case of scattering processes and statically to the decaying particle 
mass in the case of decay processes; that is 


A = 



mx, 


^rii+nf 

2-^i=ni 



V decay proc. X — {xi} 
otherwise 


(3.1) 


Processes with at least one external state Higgs accompanied by only jets and/or pho¬ 
tons do not receive any tree-level contributions even in the presence of quarks. For these 
processes we considered the contributions from the gluon and all four massless quark flavours 
in the proton and jet definitions, that are denoted by p and j. For all other processes, we 
restricted ourselves to only gluons in the external states since the loop corrections to the 
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corresponding processes with external quarks are formally NNLO. Notice that, unless the 
loop amplitudes are finite, the virtual-virtual contribution of NNLO corrections can, in 
principle, not be computed with our implementation because loop integrals are evaluated 
only up to 0{e^) in dimensional regularisation. There is however a proposal [39] for a 
formalism avoiding the computation of the 0(e) terms of one-loop amplitudes in NNLO 
computations by capturing these terms with the one-loop insertion operator instead. 

For all processes with coloured initial states, we used the MSTW 2008 PDF set [40] 
(with name ‘MSTW20081o68cl_nf4’). For each process, we indicate the maximum scale 
variation, denoted by A^, obtained by independently multiplying the scales Hr and p.fhy the 
customary factors one-half, one and two. We also show the PDF uncertainty Apop obtained 
from the corresponding MSTW error sets. Both these quantities have been computed from a 
single run with central PDF set and scale using the exact reweighting approach implemented 
by the MGSaMC module SysCalc [41]. 

We applied the cuts listed in table 3 to all processes when applicable. 


Cut 

Constraint 

Comment 

Pt,3 

> 20 GeV 

j = gluons and massless quarks 

Pt,a 7 Pt,l 

> 10 GeV 

a = photon, 1 = any lepton 

Vj 

< 5 

r] = rapidity 

ga,gi 

ARa,a, ARpi 

< 2.5 

> 0.4 

angular separation, AR = Acj)^ -\- Arf 


Table 3. General set of cuts applied (when applicable) to all processes of tables 4-8. 

In addition to these cuts, some processes are subject to more specific changes in the simu¬ 
lation setup: 

• gg ^ ZZ^ gg —)• W^W~ (b.6, b.l2) : 

Cut on heavy boson transverse momenta: pty > 1 GeV. 

The fermionic 4-point loop integral with massive identical final state vectors V features 
an integrable singularity at pt y^O that we regulated with a technical cut at 1 GeV, 
alike what is done in the latest version of the code MGFM [42]. 

• gg^ ZW^W- (c.lO): 

Complex mass scheme with Ttop = 1-49 GeV. 

Some of the diagrams contributing to this process can have up to four loop propagators 
onshell^®, and some of the resulting four-point scalar loop integrals are unstable for 
kinematic configurations with a total invariant mass slightly above twice the top mass. 

To avoid this issue, we considered here a non-zero top quark width within the complex 
mass scheme [43, 44], but keeping the weak boson widths set to zero since they are 
onshell in the external states. Such a partial assignment of the widths within the 

example of which is the pentagon with three top quark loop propagators and one bottom quark 
one, from which the Z-boson is emitted. 
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u/c djsjh t 



Figure 1. Representative diagram for the pair-production of same-sign top quarks in the SM. 

complex mass scheme is consistent at leading-order when there is no intermediate 
resonant weak bosons. In this context, the only difference w.r.t. the narrow width 
approximation is that the top quark Yukawa coupling as well as the mass in the 
numerator of top quark propagators are complex. 

• e+e- ^ HH (e.2): 

The widths of the massive bosons are set to zero in the loop propagators but kept 
finite (and equal to the values given in table 2) in the propagators of the tree struc¬ 
tures attached to the loops. 

Contrary to the two other leptonic scattering processes presented in this section, this 
one does not involve fermionic loops but genuine weak loops. We stress that the 
complex mass scheme would be the optimal approach for accounting for finite-widths 
effect since it yields gauge-invariant results.^^ 


• pp ^ tt (f.l): 

This process is extremely rare in the SM for three reasons: it only occurs via a box 
loop with two VF-boson propagators (see Fig. 1) and it is suppressed by the square of 
CKM matrix off-diagonal elements as well as the bottom quark mass because of the 
GIM mechanism. Observing such final states would therefore be a clear indication 
of new physics. We chose the following values for the additional parameters entering 
this computation: 

Tw = 0 and the CKM matrix set to non-unity, expressed in terms of the Wolfenstein 
parametrization [46], with A = 0.2253, A = 0.808, p = 0.132 and rj = 0.341. 

• Decay processes (g.x): 

Processes g.l, g.4 and g.7 are computed completely inclusively. 

For processes g.2, g.3, g.5 and g.7, all pairs of identical final state objects are 
restricted to have an invariant mass /aa > 10 GeV and all decay products must 
have an angular separation of ARj/^j/a > 0.4. 

^^This application of the complex mass scheme within EW loop computations will be made possible soon 
within the MGSaMC framework [45]. Notice however that, being UV-finite, the case of loop-induced 
processes is considerably simpler for what concerns the application of the complex mass scheme. 
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In tables 4-8, we denote by a dagger (^) processes whose inclusive cross-section has never 
been published and is, to the best of our knowledge, reported here for the first time. We 
prefix by a star (*) processes which are not readily available in the main public simulation 
tools MCFM [42], VBFNLO [47] or HPAIR [48], 


Process Syntax Cross section (pb) ^pdf Ref. 


Single boson + jets 






v^ = 

13 TeV 


a.l 

pp^ H 

P 

P 

> 

h 

[QCD] 

17.79 ±0.060 

+31.3% +0.5% 
-23.1% -0.9% 

[49] 

a. 2 

pp->Hj 

P 

P 

> 

h 

j [qcD] 

12.86 ±0.030 

+42.3% +0.6% 
-27.7% -0.9% 

[49] 

a. 3 

pp^Hjj 

P 

P 

> 

h 

j j QED=1 [QCD] 

6.175 ±0.020 

+61.8% +0.7% 
-35.6% -0.9% 

[49] 

*a.4 

99 ^ Zg 

g 

g 

> 

Z 

g [QCD] 

43.05 ± 0.060 

+43.7% +0.7% 
-28.4% -1.0% 

[34] 

*a.5 

gg^Zgg 

g 

g 

> 

Z 

g g [QCD] 

20.85 ±0.030 

+64.5% +1.0% 
-36.5% -1.1% 

[50] 

ta.6 

99^19 

g 

g 

> 

a 

g [qcD] 

75.61 ±0.200 

+73.8% +0.7% 
-41.6% -1.1% 

[-] 

Ia.7 

99 199 

g 

g 

> 

a 

g g [QCD] 

14.50 ±0.030 

+76.2% +0.6% 
-40.7% -1.0% 

[-] 


Table 4. Inclusive cross-sections for loop-induced single electroweak boson production in asso¬ 
ciation with up to two jets/gluons. A star (*) prefixes processes not readily available in the tools 
MCFM, VBFNLO or HPAIR. A dagger (I) prefixes processes whose inclusive cross-section is 
reported here for the first time. See text for details. 


Process 

Syntax 



Cross section (pb) 

A/i ApDp 

Ref. 

Double bosons ± jet 







= 13 TeV 


b.l 

pp^ HH 

P 

P 

> 

h 

h 

[qcD] 

1.641 ± 0.002 • 10-2 

+30.2% +1.1% 
-21.7% -1.2% 

[48] 

b.2 

pp^ HHj 

P 

P 

> 

h 

h 

j [qcD] 

1.758 ± 0.003 • 10-2 

+45.7% +1.2% 
-29.2% -1.2% 

[51] 

*b.3 

pp^H'^j 

P 

P 

> 

h 

a 

j [qcD] 

4.225 ± 0.006 • 10-3 

+38.6% +0.4% 
-25.9% -0.7% 

[52] 

*b.4 

gg^HZ 

g 

g 

> 

h 

z 

[qcD] 

6.537 ± 0.030 • 10-2 

+29.4% +1.0% 
-21.3% -1.1% 

[53] 

*b.5 

gg^HZg 

g 

g 

> 

h 

z 

g [qcD] 

5.465 ± 0.020 • 10-2 

+46.0% +1.2% 
-29.4% -1.3% 

[52] 

b.6 

gg^ZZ 

g 

g 

> 

Z 

z 

[qcD] 

1.313 ±0.004 

+27.1% +0.7% 
-20.1% -1.0% 

[42] 

*b.7 

99 ZZg 

g 

g 

> 

Z 

z 

g [qcD] 

0.6361 ±0.002 

+45.4% +1.0% 
-29.1% -1.2% 

[54] 

b.8 

99^ Z-f 

g 

g 

> 

Z 

a 

[qcD] 

1.265 ± 0.0007 

+30.2% +0.6% 
-22.2% -1.0% 

[42] 

*b.9 

99 Z')g 

g 

g 

> 

Z 

a 

g [qcD] 

0.4604 ± 0.001 

+43.7% +0.8% 
-28.4% -1.1% 

[55] 

b.lO 

99-^11 

g 

g 

> 

a 

a 

[qcD] 

5.182 ±0.010 • 10+2 

+72.3% +1.0% 
-43.4% -1.3% 

[42] 

*b.ll 

99 119 

g 

g 

> 

a 

a 

g [qcD] 

19.22 ±0.030 

+59.7% +0.7% 
-35.7% -1.0% 

[56] 

b.l2 

gg^W+W- 

g 

g 

> 

W4 

• w- [qcD] 

4.099 ±0.010 

+26.5% +0.7% 
-19.7% -1.0% 

[57] 

*b.l3 

gg^W+W-g 

g 

g 

> 

W4 

• w- g [qcD] 

1.837 ±0.004 

+45.2% +0.9% 
-29.0% -1.1% 

[58] 


Table 5. Inclusive cross-sections for loop-induced double electroweak boson production in asso¬ 
ciation with up to one jet/gluon. A star (*) prefixes processes not readily available in the tools 
MCFM, VBFNLO or HPAIR. A dagger (1) prefixes processes whose inclusive cross-section is 
reported here for the first time. See text for details. 
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Process 

Syntax 

Cross section (pb) 

A/1 ApDp 

Ref. 

Triple bosons 



= 13 TeV 


*c.l 

pp-^HHH 

p p > h h h [QCD] 

3.968 ±0.010 • 

10-5 

-1-31.8% -1-1.4% 
-22.6% -1.4% 

[59] 

'l'c.2 

gg^HHZ 

g g > h h z [QCD] 

5.260 ± 0.009 • 

10-5 

4-31.2% 4-1.3% 
-22.2% -1.3% 

[-1 

R.S 

99 HZZ 

g g > h z z [QCD] 

1.144 ± 0.004 • 

10-* 

4-31.1% 4-1.2% 
-22.2% -1.3% 

[-1 

R.d 

gg^HZj 

g g > h z a [QCD] 

6.190 ± 0.020 • 

10-5 

4-29.3% 4-1.0% 
-21.2% -1.2% 

[-1 

R.S 

pp^ 

p p > h a a [QCD] 

6.058 ± 0.004 • 

10-5 

4-30.3% 4-1.1% 
-21.8% -1.3% 

[-1 

*c.6 

gg^HW+W- 

g g > h w+ w- [QCD] 

2.670 ± 0.007 • 

10-4 

4-31.0% 4-1.2% 
-22.2% -1.3% 

[60] 

^c.7 

99 -t ZZZ 

g g > z z z [QCD] 

6.964 ± 0.009 • 

10-5 

4-30.9% 4-1.2% 
-22.1% -1.3% 

[-1 

tc.8 

99 ZZ-f 

g g > z z a [QCD] 

3.454 ±0.010 • 

10-5 

4-28.7% 4-0.9% 
-20.9% -1.1% 

[-1 

*c.9 

gg^Zj-f 

g g > z a a [QCD] 

3.079 ± 0.005 • 

10-4 

4-28.0% 4-0.7% 
-20.9% -1.0% 

[61] 

tc.lO 

gg^ZW+W- 

g g > z w+ w- [QCD] 

8.595 ± 0.020 • 

10-5 

4-26.9% 4-0.6% 
-19.5% -0.6% 

[-1 

tc.l2 

99 ^W^W~ 

g g > a w+ w- [QCD] 

1.822 ± 0.005 • 

10-2 

4-28.7% 4-0.9% 
-20.9% -1.1% 

[-1 

Table 

6. Inclusive cross 

-sections for loop-induced triple electroweak boson production. 

A star 

(*) prefixes processes not readily available in the tools MCFM, VBFNLO 

or HPAIR. A dagger 

(^) prefixes processes whose inclusive cross-section is reported here for the first time. See text for 

details. 







Process 

Syntax 

Cross section (pb) 

Ap ApDp 

Ref. 

Selected 2 — )• 4 


= 13 TeV 


td.l 

pp Hjjj 

p p > h j j j QED=1 [QCD] 

2.519 ±0.005 


4-75.1% 4-0.6% 
-39.8% -0.6% 

[62] 

*d.2 

pp HHjj 

p p > h h j j QED=1 [QCD] 

1.085 ±0.002 

• 10-2 

+62.1% +1.2% 
-35.8% -1.3% 

[63] 

td.S 

pp — > HHHj 

p p > h h h j [QCD] 

4.981 ±0.008 

• 10-5 

+46.3% +1.4% 
-29.6% -1.4% 

[-1 

td.3 

pp HHHH 

p p > h h h h [QCD] 

1.080 ±0.003 

•10-" 

+33.3% +1.7% 
-23.4% -1.7% 

[-1 

d.4 

gg -s> e+e-g+g- 

g g > e+ e- mu+ mu- [QCD] 

2.022 ±0.003 

• 10-5 

+26.4% +0.7% 
-19.4% -1.1% 

[64] 

td.fi 

pp HZjj 

g g > h z a g [QCD] 

4.950 ± 0.008 

• 10-5 

+45.8% +1.2% 
-29.3% -1.3% 

[-1 

Non-hadronic processes 


= 500 GeV 

no PDF 


*e.l 

e+e ggg 

e+ e- > g g g [QED] 

2.526 ±0.004 

• 10-5 

+31.2% 

-22.0% 

[65] 

te.2 

e+e- HH 

e+ e- > h h [QED] 

1.567 ±0.003 

• 10-5 


[-1 

'I’e.S 

e+e“ — >■ HHgg 

e+e- > h h g g [QED] 

6.629 ±0.010 

• 10-14 

+ 19.2% 

-14.8% 

[-1 

*e.4 

77 — > 

a a > h h [QED] 

3.198 ±0.005 

• 10-4 


[66] 

Miscellaneous 



= 13 TeV 


I'Ll 

pp —)■ tt 

p p > t t [QED] 

4.045 ± 0.007 

• 10-45 

+0.2% +0.9% 
-0.8% -1.0% 

[-1 


Table 7. Inclusive cross-sections for various 2 —> 4 processes as well as processes with non-hadronic 
initial states. A star (*) prefixes processes not readily available in the tools MCFM, VBFNLO or 
HPAIR. A dagger (^) prefixes processes whose inclusive cross-section is reported here for the first 
time. See text for details. 
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Process Syntax Partial width (GeV) Ref. 

Bosonic decays 


g-1 

H^jj 

h 

> 

j 

j 

[qcD] 

1.740 ± 0.0006 • 10-* 

[49] 

*g.2 

H ^jjj 

h 

> 

j 

j 

j [QCD] 

3.413 ±0.010 

• 10-* 

[49] 

^g.3 

H ^jjjj 

h 

> 

j 

j 

j j QED=1 [QCD] 

1.654 ±0.004 

• 10-* 

[-] 

g.4 

i7 —> 77 

h 

> 

a 

a 

[QED] 

9.882 ±0.002 

• 10-® 

[67] 

^g.5 

H -)■ jj 

h 

> 

a 

a 

j j [QCD] 

7.448 ± 0.030 

• 10-13 

[-] 

^g-7 

H —>■ 7777 

h 

> 

a 

a 

a a [QED] 

1.546 ±0.006 

• 10-14 

[-] 

*g.8 

Z^ggg 

Z 

> 

g 

g 

g [QCD] 

3.986 ±0.010 

• 10-® 

[34] 

Table 8. 

Partial decay widths for selected decay processes occurring via 

loops only. A 


(*) prefixes processes not readily available in the tools MCFM, VBFNLO or HPAIR. A dagger 
(^) prefixes processes whose inclusive cross-section is reported here for the first time. See text for 
details. 


4 A close-up on Higgs production 

The purpose of this section is to illustrate the capabilities of MGSaMC in the context of 
loop-induced processes and to present additional validation material. All results in this 
section are for the LHC 13 TeV using generation level cuts and SM parameters identical to 
those of sect. 3. We focus here on Higgs production via gluon fusion and compare it to the 
Effective Field Theory (EFT) approach. This theory approximates the loop by the effective 
operator: 

= (4.1) 

where H is the Higgs field, is the gluon field strength tensor and C is the Wilson 
coefficient (known up to N^LO) [68-70]. This approximation was recently used in the 
context of the computation of the Higgs cross-section at N^LO QCD accuracy, exhibiting 
a scale uncertainty of 2% only [6]. At this level of precision it is important to consider the 
corrections to the EFT approach, by including effects due to the finite top quark mass [71- 
73] as well as the interference between top and bottom quark loops [74, 75]. 

The importance of these contributions is already significant at leading order accuracy. 
In table 9, we compare the LO cross-section for the Effective Field Theory approach (EFT)^^ 
with the SM exact loop-induced process in two scenarios: first with massive bottom quark 
(labelled ‘LI’) and with massless bottom quark (labelled ‘NoB’) where the contribution of 
the bottom quark loop is vanishing (both when interfering with the top quark loop and 
when squared against itself). For the zero jet multiplicity, we have a perfect agreement 
between the ‘NoB’ case and the ‘EFT’ one because the effective field theory is valid when 
= mu < 2mt^^ In presence of radiation, the total energy of the event can lie outside 

use the built-in ‘heft’ model of MGSaMC [76]. 

^®Thls very good agreement is also due to the presence of mass correction factors in the ‘heft’ model. 
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process 

EFT 

Exact loop-induced (LI) 

Exact loop-induced mb = 0 (NoB) 

99 

19.996(4) pb 

17.79(6) pb 

19.94(4) pb 

pp hj 

13.41(2) pb 

12.86(4) pb 

13.24(4) pb 

pp hjj 

6.31(2) pb 

6.18(2) pb 

6.13(1) pb 


Table 9. Comparison of the cross-section for loop-induced Higgs production computed in different 
setups. The first column reports the prediction from the Effective Field Theory (EFT) limit using 
the ‘heft’ model [76]. The cross-sections of the second column (LI) are computed by directly 
integrating the loop-induced diagrams, hence keeping all mass effects. The setup of the last column 
(NoB) is identical to that of the second one, except for the fact that the contribution from bottom 
quark loops is removed. We do not include the vector boson fusion contribution in the cross-sections 
reported for the process pp —>■ Hjj. 


the validity range of the EFT theory and this explains the difference observed between these 
two schemes. Fig. 2 (left part) shows the Higgs pp distribution for various computational 
setups and it clearly demonstrates that the EFT breaks down for large Higgs transverse 
momenta. 

Another effect emphasised in table 9 is the importance of the interference of the top 
and bottom quark loops. For the zero jet multiplicity case, this interference is negative 
because of the absorptive part of the bottom quark loop which is below threshold, unlike 
the purely real top quark loop. In presence of additional jets, the interference can be either 
positive or negative, leading to a milder effect on the cross-section but affecting the shape of 
the distribution. In Fig. 2 (right part), we show the interference contribution as a function 
of the Higgs transverse momentum (for the process pp —t Hj). At low pp, we recover 
the behaviour of the 0-jet multiplicity case featuring a negative interference effect which 
contributes to 10% of the total cross-section. Conversely, at high pp, the interference of the 
bottom and top quark loop diagrams is constructive and its relative contribution tends to 
a constant equal to the ratio of the corresponding Yukawa couplings: = 2.7%. 

The left panel of fig. 2 features a fourth prediction labelled ‘RWGT’, short for reweight¬ 
ing. This prediction corresponds to the differential cross-section predicted by the Effective 
Field theory where the weight of each event has been rescaled by with Mpj being 

the loop-induced matrix element for the helicity configuration h computed for the kine¬ 
matic configuration of the event considered. With infinite statistics, the prediction based 
on reweighting must be identical to the one obtained from a direct integration of the loop- 
induced matrix element. We check that the ratio between these two predictions (shown in 
the bottom inset) is always compatible with one, given the statistical uncertainty whose 
one-sigma variation is given as the yellow band. This represents a non trivial validation of 
MadEvent phase-space integration for loop-induced processes. 

As it is the case for tree-level computations, MGSaMC can perform matched and 
merged computations in the context of loop-induced processes. The MLM merging (either 
/ct-MLM or shower kp [77]) is available and fully automatic when linked to PythiaG [78]. 
As an example of such computation, we present plots for Higgs production merged with up 
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Figure 2. Differential distribution of the Higgs transverse momentum. The left panel presents the 
prediction obtained within the three different computational setups LI, EFT and NoB introduced 
in table 9. The red curve, labelled ‘RWGT’, is obtained by reweighting events generated from the 
effective theory approach with the exact loop-induced matrix element including all top- and bottom 
quark effects (therefore strictly equivalent to the ‘LT prediction in the infinite statistics limit). The 
bottom inset of the left panel shows the ratios between the various predictions, overlaid (in yellow) 
by the one sigma statistical error band. The right panel shows the contribution of the interference 
between the top and bottom quark loops, and its bottom inset indicates its contribution relative to 
the total cross-section. 


Differential Jet Rate 1 >0 Differential Jet Rate 2 >1 



Figure 3. First two Differential Jet Rate (DJR) histograms for exact loop-induced Higgs produc¬ 
tion predictions matched to Pythia6 and merged with up to two additional jets. The bottom insets 
show the ratio of the prediction using exact loop-induced matrix elements to the one performed 
within the Effective Field Theory (EFT) framework. 


to two additional jets in the /ct-MLM scheme. In Fig. 3, we present the differential jet 
rate plot (DJR< 7 v>) which corresponds to the scale for which the kx clustering algorithm 
switches from an N to N-1 jet description, that is the scale associated to the emission. 
This observable is most sensitive to the matching/merging since it is directly related to 
the variable disentangling matrix-element emissions from parton shower ones. The parton 
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Figure 4. Higgs transverse momentum differential distribution in a matched/merged (fcr-MLM, 
Qmatch = 50 GeV) prediction with up to two jets simulated with matrix elements. We present pt 
distributions for both the Higgs and the hardest (p^-ordered) photon into which it decays. 

shower is not allowed to radiate in the phase-space region above the merging scale (set 
here to Qmatch = 50 GeV as in [79]) since this region is already popnlated by partonic 
events generated nsing matrix elements of higher mnltiplicity. Therefore, having a smooth 
transition at Qmatch for the snm of samples (distribntion in black) is an important check 
of the qnality of the matching/merging procednre. For DJRi, a small discontinnity can be 
observed, originating from the fact that the bottom and top qnark interference effect is not 
acconnted for by the parton shower. 

Finally, fig. 4 shows the normalised distribntion of the Higgs transverse momentnm 
obtained from the same matched and merged sample. This plot also presents the transverse 
momentnm of the hardest photon when considering the Higgs decaying into two photons 
(simnlated in the narrow-width approximation). Notice that the discrepancy between the 
EFT and the loop-indnced prediction is only dne to differences in the prodnction mechanism 
since the (loop-indnced) matrix element for any 2-body decay is a constant over phase-space. 

5 Loop-induced processes within BSM models: the 2HDM example 

Loop-indnced processes simnlation with MG5 aMC is not limited to the Standard Model 
bnt can also be extended to a large class of BSM models. The framework snpports mod¬ 
els following the NLO UFO format [80-82] which can be antomatically created by Feyn- 
Rnles/NLOCT [83, 84]. Since loop-indnced processes are hnite, it is not reqnired to have 
UV connterterms specihed in the NLO UFO model, therefore allowing to compnte loop- 
indnced processes even for models for which antomatic UV renormalisation is cnrrently not 
possible. The compntation of the R 2 Feynman rnles is however still necessary (see sect. 2.1). 

^“^There are two exceptions to this rule: first, for the highest multiplicity sample, the shower can radiate 
up to the scale of the softest matrix element jet and secondly a jet can be excluded from the matching if an 
emission is due to a diagram which does not have a counterpart in the parton shower description. In the 
latter case, the event is seen by the matching procedure as having one less jet, since the corresponding jet 
emission is excluded. This explains why the dashed blue curve (corresponding to the single jet multiplicity) 
does not vanish above the matching scale for the DJR 2 plot (left panel of fig. 3). 
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Type 

tan (3 

aj-K 

rriho 

rriffo 

m^o 

mH± 

m?2 

B1 

H 

1.75 

-0.1872 

125 

300 

441 

442 

38300 

B2 

I 

1.20 

-0.1760 

125 

200 

500 

500 

-60000 

B3 

H 

1.70 

-0.1757 

125 

350 

250 

350 

12000 


Table 10. Benchmark parameter sets for the two Higgs Doublet Model (2HDM) used for the 
computations presented in this section. All masses are given in GeV. The top quark mass is set to 
173 GeV and the bottom quark mass to 4.75 GeV. 



gg ^ Zh^ 

gg ^ ZH^ 

gg ^ 

B1 

B2 

B3 

1 1 q R +28.9% +1.0% 
IIO.U _2i.2% -1,2% 

or: CQ +29.9% +1.4% 
00.oy _2 i.4% _i,i% 

IfiO q +28.1% +1.4% 
19.9% o,5% 

fjon A +29.6% +1.2% 

UOZ.4 _21.5% -1,2% 

1 +30.1% +1.3% 

1040 _2i.8% _i.3% 

n SQRS +31.2% +1.5% 
0-8068 22,3% -1.6% 

n R9nq +32.5% +i.9% 

U.UZUO _23,o% -1.9% 

n 8R1 A +33.0% +2.0% 
U.5Di4 _23,3% _2.o% 

1 *^1 7 +28.4% +1.0% 
-20.8% -1.0% 


Table 11. Inclusive cross sections (in fb) for loop-induced Z boson production in association 
with each neutral Higgs of the 2HDM model at the 14 TeV LHG. We report results for the three 
benchmark parameters sets introduced in table 10. The first uncertainty (in percent) refers to scale 
variations by factors ^, 1 and 2 while the second one refers to PDF uncertainty. Both are computed 
by SysCalc at no additional computational cost. 


B1 

B2 

B3 


gg 


H+H- 


-1-34.0% -1-2.2% 
-23.8% -2.2% 
-1-34.6% -1-2.4% 
-24.1% -2.4% 

n S +32.8% +1.9% 
U.DiO 23,2% -1.9% 


.2334 

'.7011 


qq 


H+H- 


0 7669 +5-9% +1-1% 

U. (UUy 

0 4406 +1.4% 

U.'lUUU 

9 079 +4.6% +0.9% 
Z.U^Z 


Table 12. Inclusive cross sections (in fb) for the production of a pair of 2HDM charged scalars 
at the 14 TeV LHG for the three benchmark parameters sets introduced in table 10. The first 
uncertainty (in percent) refers to scale variations by factors ^,1 and 2 while the second one refers 
to PDF uncertainty. Both are computed by SysCalc at no additional computational cost. 


In this section, we present results in the two Higgs Doublet Model (2HDM) [85] for 
the three benchmark parameter sets presented in table 10 and introduced in refs. [27, 86]. 
Ref. [27] computes the cross-section of the associated production of a Z boson with all 
the 2HDM scalars using a reweighting approach. As a validation of our method, we have 
reproduced the computation presented in table 6 of ref. [27]. The only two differences 
with respect to the parameters used in sect. 3 is the c.o.m. energy set to 14 TeV and the 
renormalisation and factorisation scales set to the partonic invariant mass. Our results are 
reported in table 11 and show perfect agreement. 

In table 12, we present results for the production of a pair of charged Higgs via gluon 
fusion (see Fig. 5) and we compare it to the tree-level channel with initial-state quarks. 
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Figure 5. Representative diagrams for the production via gluon fusion of a pair of charged Higgs 
in the two Higgs doublet model. 


production of a pair of charged Higgs at 14 TeV LHC production of a pair of charged Higgs at 14 TeV LHC 




Figure 6. Differential distributions at the 14 TeV LHC for the production of a pair of charged 
Higgs. We present separately the loop-induced gluon initiated production (solid line) and the 
tree-level quark initiated one (dashed line). 


Even if the loop-induced contribution is here formally NNLO, it is numerically important. 
Moreover, the loop-induced contribution differs significantly in its kinematic distribution, as 
illustrated in Fig. 6 where we show the system invariant mass and the transverse momentum 
distribution of each partonic subprocess. For all benchmark points, the distributions are 
softer for the loop-induced production channel because of the differences in the relevant 
parton density functions. 

6 Conclusion 

We have presented a new public extension of MadGraph5_aMC@NLO capable of auto¬ 
matically simulating loop-induced processes within both the Standard Model and BSM 
theories. Our implementation does not rely on any predefined set of processes and instead 
generates on demand an optimised code tailored to the loop-induced process of interest. 
This is achieved by combining the capabilities of MadLoop for the generation of loop- 
induced matrix elements with those of MadEvent for phase-space integration. 

Particular emphasis is put on run time and we have optimised both aforementioned 
codes for the purpose of loop-induced processes whose matrix-elements are considerably 
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slower to evaluate than tree-level ones. We have first determined that, within our frame¬ 
work, the most efficient approach is to employ the OPP reduction method for a given 
helicity configuration, which we pick for each phase-space point according to a probability 
set by importance sampling. A second improvement is the computation of the loop matrix 
elements by projecting its colour structures onto the colour-flow basis. This also provides 
the associated leading colour information necessary for event generation and required by 
parton shower programs. Finally, we have restructured MadEvent code so as to allow for 
an arbitrary parallelisation of the computation. 

Thanks to these developments, it is possible to compute all 2 —)■ 2 loop-induced pro¬ 
cesses on a laptop and all 2 —)■ 3 ones (and most 2 —)■ 4) on a small size cluster. As a 
proof of this statement, we have computed a large list of cross-sections in the Standard 
Model including all loop-induced processes up to 2 —)■ 3 and a couple of 2 —)■ 4 processes as 
well as some decay and non-hadronic scattering processes. The implementation was care¬ 
fully validated against results in the literature for both integrated distributions and local 
phase-space points. 

Our approach is fully differential and compatible with standard matching and merging 
procedures like the MLM approach supported by MGSaMC; this is emphasised by phe¬ 
nomenological results and distributions presented for the case of Higgs production. Finally, 
the flexibility and generality of MGSaMC regarding the support of BSM models is main¬ 
tained and we illustrated this with results for various loop-induced processes within the two 
Higgs doublet model. 

This work fills a gap in the spectrum of modern tools capabilities and it opens the 
door to many applications whose exploration will no longer be hindered by the technical 
difficulties of simulating loop-induced processes. 
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A Manual 


In this appendix, we present minimal information on how to generate events for a given 
loop-induced process. We then proceed with the description of the new options introduced 
in MG5 aMC version 2.3.0 and their particular application in the context of loop-induced, 
tree-level and NLO computations. 

A.l Basic commands 

The steps needed to generate loop-induced processes are very similar to the ones needed 
to generate a tree-level and/or an NLO computation within the MGSaMC framework. For 
the reader familiar with that environment, it should be enough to stress that for processes 
without any Born diagrams, the ‘generate’ command suffixed with the NLO syntax (he. 
adding ‘ [QCD] ’ at the end of the process definition) will automatically consider the corre¬ 
sponding loop-induced process. From now on, we adopt a pedantic approach and do not 
assume any prior knowledge of the tool. 

The MGSaMG framework is a meta-code generating optimised numerical programs for 
the simulation of user-defined scattering/decay processes [18]. The code offers an interactive 
interface, that can be started with the script <MG_install_path>/bin/mg5_aMC, where the 
user can enter commands to specify the processes/runs he is interested in.^® All commands 
and options are documented within the code, and detailed information can be obtained 
by typing ‘help COMMAND’. The list of all commands available are listed by typing ‘help’. 
Finally, a built-in tutorial can be started using the command ‘tutorial’. We now proceed to 
describing the succession of the four commands characteristic of the steering of a simulation. 

• import model <MODEL_NAME>: 

This command allows to import a new model. In the context of loop induced pro¬ 
cesses, only NLO UFO models are supported. Such models can be generated via 
FeynRules [83] thanks to the NLOCT package [84]. If this command is skipped, a 
simplified SM model (with diagonal-CKM and massless quark/leptons up to second 
generation) is loaded by default, 
examples: 

import model loop_qcd_qed_sm 
import model 2HDM_NLO 

• generate <PROCESS_DEEINITION>: 

This line corresponds to the definition of the process of interest and to the generation 
of the corresponding Feynman diagram. The process definition consists in a string 
with the initial state and final state particles (separated by ‘>’). To specify a loop- 
induced process, one must suffix the process definition with the tag ‘ [QCD] ’ (same 
as for NLO QCD computations). In the presence of this tag, MadLoop first checks 

^®A11 those steps can also be done without using the interactive interface (i.e. scripting). For more 
information on this please read: https://answers.launchpad.net/mg5amcnlo/+faq/2186 
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whether it receives tree level contributions; if it does, then a standard NLO compu¬ 
tation is performed. It is possible to use the tag ‘ [noborn=QCD] forcing MadLoop 
to build matrix elements using the square of loop diagrams (be. loop-induced mode) 
without checking for the existence of tree-level contributions, 
examples: 

generate p p > h j [QCD] 

generate p p > t t [QED] 

If one is interested in only the loop matrix element evaluation for a given kinematic 
configuration (be. not including the integration and event generation related code), 
the tag [virt=QCD] should be used instead. The subsequent ‘launch’ command then 
asks for the specification of the phase-space point and model parameters and returns 
the corresponding numerical value of the loop matrix-element. A dynamic library of 
the fortran code produced by MadLoop for the matrix element computation can then 
easily be produced by running the shell command ‘make OLP’ in the ^SubProcesses' 
directory of the process output folder, 
example: 

generate h > g g [virt=QCD] 

• output <PATH>: 

This starts the generation, in the specified path, of the numerical code for the process 

of interest. 

example: 

^^ut MY_NEW_RUN 

• launch <PATH>: 

This command starts the actual computation of the cross-section and the generation of 
events. The code first asks two preliminary questions. First, the list of external tools 
(Pythia6 [78], Delphes [87], MadSpin [88], ...) that one wishes to use in this simula¬ 
tion can be turned on. Then, the user is given the option to edit various configuration 
files, including the model parameters. This part is identical to the case of leading or¬ 
der computations and we refer to [18] for more details in this subject. Once the 
simulation is completed, various outputs are written in the folder ^<PATH>/Events' 
and a convenient html summary is generated in ^<PATH>/crossx.html'. 

example: 

launch MY_ NE W_ R UN 

A.2 Description of new MadGraph5_aMC@NLO options 

In this new version of MG5aMC (version 2.3.0), we have introduced a couple of additional 
options related to the simulation of loop-induced processes. The following options can be 
modified in the file ^input/mg5_configuration.txt' or directly via the interface using the 
‘set’ command. 
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• cluster_size [default—100]: This parameter allows to arbitrarily scale the requested 
number of jobs to be proportional to the size of the cluster. By increasing the cluster 
size, MG5aMC splits iterations into more jobs, each probing less phase-space points. 
As of now, this parameter applies for the simulation of loop-induced processes only. 

• cluster_local_path [default—None]: This parameter avoids either to transfer PDF 
sets to the cluster nodes or to read them directly on a central disk. This path should 
point to a (node specific) directory containing the associated PDF sets (either those 
from LHAPDF or built-in ones). A typical usage is to set this path to a local directory 
mirrored via cvmfs. 

• max_npoint_f or_channel [default—3]: This parameter controls what topologies enter 
the multi-channeling to be used for integrating loop-induced processes. For instance, 
when set to 4, all loop diagrams with 4 or less loop propagators seed their own channel 
of integration. This means that tree topologies obtained by shrinking up to box loop 
diagrams are considered for the multi-channeling and therefore integrated separately. 
In general, we do not observe any significant gain (when not detrimental) when setting 
this parameter to 4 or larger (except in the case of gg —)■ zz). We stress here that, 
for this parameter to take effect, it must to be modified prior the generation of the 
source code of the process considered. 

• loop_colour_f lows [default—dynamical]: The computation of partial colour sub¬ 
amplitudes (z.e. amplitudes for fixed colour flows) is turned off by default for the 
case of NLO virtual matrix elements and turned on for loop-induced matrix el¬ 
ements.. This is because in the former case it comes at the price of giving up 
loop reduction at the squared amplitude level, hence slowing down MadLoop exe¬ 
cution speed since the number of OPP reductions is no longer independent of the 
number of contributing helicity configurations. This option can however be turned 
on (before MadLoop writes out the source code for the process) since colour sub¬ 
amplitudes can be necessary for certain applications, such as NLO event genera¬ 
tion within the context of a controlled colour expansion [89] and/or Monte-Carlo 
over colours. Partial colour amplitudes (called JAMP in the code) can then be ac¬ 
cesses and combined as needed by modifying the user-defined subroutine ‘ *_ COM- 
PUTE_COLOR_FLOWS_DERIVED_QUANTITIES' present in the source code 
file ^compute_color_flows.f'. 

• loop_optimized_output [default—True]: This option corresponds to using the poly¬ 
nomial decomposition of the integrand numerator (Eq. 2.6), introduced in [13], to 
optimise the use of OPP reduction and allow to interface TIR tools. For debug¬ 
ging and validation purposes it is useful to be able to turn this option off (before 
generating the loop matrix element code) and to force MadLoop to recompute the 
complete integrand numerator for each new value of the loop-momentum specified by 
the OPP reduction procedure. This provides a strong check on the correctness of the 
computation of the polynomial coefficients. 
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value 


-1 


0 

1 


dynainical_scale_choice meaning 
default case: 

LO code: transverse mass of the 2—7-2 system resulting 
of a kx clustering 

NLO code: sum of the transverse mass divide by 2 
user defined scale, specified in the file ’setscales.f ’ 
total transverse energy of the event 




2^i=l 


Ei-pT,i 




2 

3 

4 


sum of the transverse mass 

sum of the transverse mass divide by 2 

partonic energy 


i:f.i sj’^i+pir 

5 Eil \/+ pt.i- 


' S. 


Table 13. Values supported in ^run_card.dat'’ for the new entry ‘dynamical_scale_choice’ and 
their corresponding functional form of the factorisation and renormalisation scale fip and fj^n. This 
parameter applies for both leading order and next to leading order computations. 


On top of the parameters above which control the way the code is generated and then 
handled by the cluster, we also introduced new entries in ^run_ card. dat\ which can be 
accessed and modified at any time, including after code generation. 

• nhel: This parameter was already defined in previous versions of the code but its effect 
changed in MGSaMC v2.3. ‘nhel’ can now only be set to ‘0’ or ‘1’, in which case all, 
respectively exactly one, helicity configuration(s) are/is considered for each phase- 
space point. The difference with previous versions is that the helicity configuration 
picked for each phase-space point is chosen according to a dynamical importance 
sampling Monte-Carlo method. The default value for this parameter is taken to be 
‘1’ for loop-induced process and ‘0’ for tree level computations (this parameter is not 
available for NLO computations). We stress that when setting this parameter to ‘0’ for 
loop-induced processes, it is more efficient to select the tensor integral method PJFry 
or IREGI as the preferred reduction by modifying the parameter ‘MLReductionLib’ 
of the ^MadLoop_card.daV configuration file. 

• dynamical_scale_choice: When a dynamical scale is chosen, this parameter selects 
its functional form (c./. table 13 for a list of predefined functional forms). 

• survey_splitting: For the phase-space integration of loop-induced processes, this 
parameter allows to define the number of nodes assigned to the integration of one 
given channel integration during the survey/gridpack generation. Notice that the 
parallelisation of the second stage of the simulation (which completes event generation 
and is referred to as the refine step) is controlled by the parameter cluster_size 
instead. 


- 25 - 













• job_strategy: This parameter controls the parallelisation strategy for the simulation 
of high multiplicity tree level simulations (especially important for an efficient han¬ 
dling of the simulation of the highest multiplicity sample in the context of an MLM 
merging gridpack generation.). The three possible integer values taken by this option 
are: 


— ‘O’: The original behaviour already adopted in previous versions where each sub¬ 
mitted jobs performs successively the integration of two channels of integration. 

— ‘1’: This mode changes the behaviour for the highest multiplicity sample, where 
the associated jobs run a single channel of integration. The number of jobs 
submitted in this case is therefore increased by two, hence reducing the running 
time for this sample by the same factor. 

— ‘2’: This mode uses the same algorithm as the one introduced for loop-induced 
computations for the highest multiplicity sample and uses the mode ‘1’ for the 
next-to-highest multiplicity sample. 

The default value of some of the Wun_card' parameters are now dynamically chosen 
depending on the process considered, so as to better reflect what is typically expected in 
this case. These parameters remain accessible and can be modified by the user; only their 
default value is changed. The parameters for which the default value can differ depending 
on the process considered are: 

• energy of the beam: for electron-positron collision, the default energy is set to 1 
TeV (500GeV per beam). The default collision energy for all other processes is 13 
TeV. 

• PDF type: if the particles in the initial states of the process are constituted of quarks, 
gluons or photons then the proton PDFs are used. Otherwise, this parameter is set 
to the fixed energy scheme without PDF. 

• MonteCarlo over helicity: Monte-Carlo over helicity is turned off for tree-level 
simulations and turned on for loop-induced ones. This options is not available for 
NLO computations. 

• maxjetflavour/asrwgtflavour: set accordingly to the number of quark flavours 
appearing in the initial states. 

• cuts: for (LO) width computation, all cuts are turned off by default. Unlike previous 
versions of the code, the cuts defined in the run_card.dat are now applied to the 
computation of the partial width as well. 

• matching parameter: At LO, if all the processes considered only differ by their jet 
multiplicity, the MLM matching/merging (ickkw=l) scheme is turned on by default 
with a parton level cut (‘xqcut’) of 30 GeV. 
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Finally, a new MadSpin option allows to simulate arbitrary decays (including three- 
body and loop-induced decays) but without any spin correlation (between production and 
decay) and without BreitWigner smearing. The detailed description of this new functional¬ 
ity is deferred to a future work, and here we limit ourselves to mentioning that this feature 
can be used by adding the line ‘set spinmode none’ at the beginning of the MadSpin 
configuration file. 

B Benchmark results for various loop-induced Higgs production pro¬ 
cesses 

We start here by presenting the numerical result for the evaluation of the matrix element 
for the loop-induced process gg —)• hhgg, chosen both for its complexity and importance 
as a background to double Higgs production in vector boson fusion. The SM parameters 
used in this computation are those specified in table 2 with ag set to 0.118. The kinematic 
configuration considered is: 


[GeF] 

E 

Px 

Pv 

Pz 

Psi 

= ( 500 

, 0 

, 0 

, 500 

P92 

= ( 500 

, 0 

, 0 

, -500 

Phi 

= ( 148.556611322403 

, -20.0350647739655 

, 36.3342710976818 

, -68.7203295313605 

Phi 

= ( 322.824972379014 

, -94.1433807610584 

, -273.713728753781 

, 69.3453772222178 

P 95 

= ( 138.118065647892 

, -95.9848700486689 

, -88.5772755257372 

, 44.9173801606153 

P96 

= ( 390.500350650690 

, 210.163315583693 

, 325.956733181836 

, -45.5424278514726 


We report here the squared loop-induced matrix element^® computed for the phase-space 
point above, summed over all helicity and colour configurations. 


gg -1 hhgg 

[GeV-^] 

\M^LI)^2 

1.3033633142042775e-12 


This particular computation is performed in quadruple precision and stability tests show 
that all 17 double precision digits are numerically stable. As for NLO computations, the 
standalone MadLoop output can be used for crosschecking an independent implementation 
of the calculation (see how in sect. A.l). 

We now turn to listing MadLoop performances in table 14. On top of the process 
gg —)• hggg, we also consider gg —>■ hh, gg —)• hhg and gg —)• hggg so as to reflect the scaling 
of MadLoop timings with the multiplicity and number of Feynman diagrams. The profiling 
presented for these processes, as well as for any other, can be automatically reproduced 
by running the command ‘check profile <process_def inition>’ from the MG5aMC in¬ 
teractive interface. The timings indicated do not include any numerical stability test nor 
do they account for the fraction of points for which it is necessary to use quadruple preci¬ 
sion arithmetics. The percentages in parenthesis specify the fraction of the running time 

^®Notice that the customary one-loop matrix element prefactor (^) irrelevant in the context 

of a finite loop matrix element. The dependence on the renormalisation scheme/scale only comes through 
the running of the Os. 
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gg -i- hh 

gg hhg 

gg hhgg 

gg -t hggg 

loop Feynman diag. 

16 

108 

952 

2040 

^ topologies 

8 

54 

380 

540 

helicity conhg. 

2 

8 

16 

32 

Generation time 

8.7s 

21s 

269s 

lh36m 

Output code size 

0.5 Mb 

0.7 Mb 

1.8 Mb 

3.2 Mb 

Runtime RAM usage 

4.7 Mb 

20.5 Mb 

102 Mb 

240 Mb 


Timing for the computation of a single helicity configuration 


OPP as in CutTools 

TIR as in IREGI 

TIR as in PJFry 

TIR as in Golem95 

2.6ms (19%) 
17.5ms (3%) 
3.2ms (15%) 
15.1ms (3%) 

40.7ms (16%) 
1.14s (0.6%) 
190ms (4%) 
615ms (1.2%) 

795ms (13%) 
65s (0.17%) 
28s (0.38%) 
16s (0.67%) 

1.03s (15%) 
79s (0.18%) 
29s (0.50%) 
19s (0.75%) 

Timing for the computation summing all helicity conhgurations 

OPP as in CutTools 

TIR as in IREGI 

TIR as in PJFry 

5.2ms (18%) 
18.4ms (5%) 
3.8ms (25%) 

328ms (15%) 
1.19s (4%) 
243ms (21%) 

14.7s (19%) 
68.2s (2.6%) 
33s (6%) 

33s (14%) 
82.0s (6%) 
38.2s (14%) 


Table 14. Performances of MadLoopS for the computation of various benchmark processes, tested 
with the gfortran compiler, no optimizations, v4.8.2, on an i7, 2.7 GHz CPU. The percentages 
in parenthesis correspond to the fraction of time spent for the computation of the polynomial 
coefficients (the complementary time is entirely spent in the loop reduction). The number of 
helicity configurations reported includes only those whose contribution is not analytically zero. 

spent for the computation of the polynomial coefficients of the numerator integrand. The 
complementary time is entirely spent in the loop reduction algorithm. Table 14 illustrates 
various points discussed in the technical sect. 2.1, especially the fact that only the OPP loop 
reduction time scales with the number of helicity computed. We refrained from showing 
the timing using Golem95 reduction for the case where the matrix element is summed over 
helicity conhgurations, because the recycling of tensorial coefficients is not yet implemented 
for this tool. In general, even though TIR reduction allows for faster timing when summing 
over helicity conhgurations, it remains much larger than for the computation of a single 
helicity conhguration, hence our default approach of using OPP reduction in conjunction 
with a Monte-Carlo over helicity conhgurations. 

The number of topologies refers to the number of different sets of loop propagator 
denominators present in the computation. In the computation of the virtual matrix element 
for NLO predictions, the integrand numerators of all loops sharing the same denominator 
topology can be added together before being reduced, hence greatly diminishing the number 
of necessary OPP reductions. This grouping of topologies is not applicable for loop-induced 
computation when using OPP reduction, but TIR can beneht from it since tensor integrals 
can be recycled across loops sharing the same denominator topology. 
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Reduction tool 

Max. nioop prop. 

Max. rank 

Complex masses 

GutTools 

lot 

^loop prop. 1 

yes 

IREGI 

7t 

7t 

yes 

PJFry 

5* 

^loop prop 

no 

Golem95 

6 

max(6, nioop prop. + 1) 

yes 


This limitation is not intrinsic to the reduction tool, and is only parametrical 
so that it is trivial to increase if proven necessary. 

Pentagons in PJFry are formally supported but typically too unstable for integration. 


Table 15. Limitations of the different reduction methods interfaced to MadLoop. 

It is interesting to note that the rednction time in TIR is almost the same between the 
processes gg —)• hhgg and gg —)• hggg, even thongh the latter has twice as many diagrams. 
This is mainly becanse the rednction time is dominated by the rednction of topologies with 
maximal tensorial rank {rmax = 6 here), of which there is the same nnmber in these two 
processes. 

Generation time, ontpnt code size and RAM nsage show that MadLoop is light-weight 
and practical for the compntation of loop-indnced amplitndes with np to at least 5 external 
legs, and well-snited for cross-checking other codes for more complicated processes. 

The varions rednction methods compared in table 14 have different ranges of applica¬ 
bility, which we snmmarise in table 15. MadLoop checks these constraints at rnntime for 
each loop and for each available rednction tool before nsing the first one applicable. It is 
possible that a mixtnre of different rednction methods is nsed within a single loop matrix 
element compntation. Notice that, becanse of a certain class of hexagons which are not 
snpported, PJFry is by defanlt limited to rank-5 loops. We have disabled that limitation 
for the present benchmark since the timing is unaffected. Also, IREGI is the only rednction 
tool capable of redncing loops with tensorial rank exceeding the nnmber of loop propagators 
by more than one nnit. 
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